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ABSTRACT 

Context. RS Ophiuchi is a recurrent nova with a period of about 22 years, timescale, consisting of a wind accreting binary system 
with a white dwarf (WD) very close to the Chandrasekhar limit and a red giant star (RG). The system is considered a prime candidate 
to evolve into an SNIa. For its most recent outbursts in 1985 and 2006, exquisite multi wavelength observational data are available. 
Aims. Further physical insight is needed regarding the inter-outburst accretion phase and the dynamical effects of the subsequent nova 
explosion in order to improve the interpretation of the observed data and to shed light on whether the system is an SNIa progenitor. 
Methods. We present a three dimensional hydrodynamic simulation of the quiescent accretion and the subsequent explosive phase. 
Results. The computed circumstellar mass distribution in the quiescent phase is highly structured with a mass enhancement in the 
orbital plane of about a factor of 2 as compared to the poleward directions. The simulated nova remnant evolves aspherically, propagat- 
ing faster toward the poles. The shock velocities derived from the simulations are in agreement with those derived from observations. 
For v RG = 20 km/s and for nearly isothermal flows, we derive a mass transfer rate to the WD of 10% of the mass loss of the RG. For 
an RG mass loss of 10 _7 M o yr _1 , we found the orbit of the system to decay by 3% per million years. With the derived mass transfer 
rate, multi-cycle nova models provide a qualitatively correct recurrence time, amplitude, and fastness of the nova. 
Conclusions. Our 3D hydrodynamic simulations provide, along with the observations and nova models, the third ingredient for a 
deeper understanding of the recurrent novae of the RS Oph type. In combination with recent multi-cycle nova models, our results 
suggests that the WD in RS Oph will increase in mass. Several speculative outcomes then seem plausible. The WD may reach the 
Chandrasekhar limit and explode as an SN la. Alternatively, the mass loss of the RG could result in a smaller Roch volume, a common 
envelope phase, and a narrow WD + WD system. Angular momentum loss due to graviational wave emission could trigger the merger 
of the two WDs and - perhaps - an SN la via the double degenerate scenario. 

Key words, stars: individual: RS Oph - stars: novae, cataclysmic variables - supernovae: general - accretion, accretion disks - 
Gravitational waves - methods: numerical 



1 . Introduction RS Oph has an orbital period of 455 days dFekel et al.ll2000h 

rN RG and WD masses of 2.3 M© and close to 1.4 M© (Anupama 

Jh ; Type la supernovae (SNe la) are cornerstones of modern cos- I2002I: IWorters et^D l2007h . respectively, and a separation be- 

mology, as a measure for cosmological distances, and they tween me components of a = 2 .68 • 10 13 cm. Th e RG mass 

are crucial building blocks of the universe, as production sites ^ however, still uncertain, much smaller values (iFekel et al.1 

of a large part of iron group elements. The recurrent nova ^ Q) have also been suggest ed. We supplement the system pa- 

RSOphiuchi (RS Oph) ( iHachisuetalJ |1999| ; | Hachisu & Kato| rameters by assuming for the RG a terminal wind velocity of 

£00iD is a candidate for their still mysterious progenitors. It Vrq = 2Q km/s in me rest frame of me RG and a mass loss rate of 

is a symbiotic-like binary star system consisting of a red gi- 1Q -7 Mp)Vr -i . Note that mass loss es of red giants are still not well 

ant (RG) and a white dwarf (WD) that undergoes a nova known. While Se aquist & Tavl or (1990) find a massloss of order 

outburst about every 22 years (|Anupama| |2002|; | Fekeletal. | lo^yr" 1 for me majority of RG in symbiotic systems, the 



|2000l iDobrzycka et al. 1996). For the most recent outbursts in scatter is significant and values of less than ip-J ^yr- 1 can be 
1985 and 2006, exquisite panchromatic observational data are found For the range of spectral types suggested Jw ort ers et al.1 



i ■ — it 1 -| | 1| 1 | 1 luunu. iui liic icuigt; ui apc^uai lypt^ suggcsicu v vvuiicis ci ai. 

available (|Shoreet al. |1996|; |Bode et al. 2006; |Sokoloski et alJ [2007) for RGs in the RS Oph system, its radius is always smaller 

2006; Das et al. 2006; O'Brien et al. 2006; Evans et al. 2006; than its Roche lobe and accret ion by the WD occurs only from 

Worters et al. 2007; Hachisu et al. 2007; Bode et al. 2007). Here me RG wind 

we present first 3D hydrodynamical simulations of RS Oph, 

from pre-outburst accretion through nova explosion, with the 

goal of improving the physical insight into the system and the 2. 3D simulations Of RS Oph 

interpretation of the unique observational data. . ~ ... , ... 

r n 2.1. Computational method 

Send offprint requests to: Rolf Walder Our n umerical simulati o n was performed using the A-MAZE- 

* In the frame of the computing project 'Cosmic Engines in Galaxies' code (Walder & Folini 2000), a parallel, block structured, 
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Fig. 1. The grid structure of the nine levels of refinement, shown 
for the entire computational domain (left) and around the accret- 
ing WD (right). While the coarser grids of level 1 through 6 are 
fixed in space, the meshes of level 7 through 9 follow the WD. 
The RG is shown in red, the accreting WD in blue. 



adaptive mesh refinement (AMR) hydrodynamical code us- 
ing Cartesian meshes and a multidimensional high-resolution 
finite- volume integration scheme . The code ha s been used for 
accretion studies before ( Zarinelli et al] 1 19951: IWalderl 119971: 
Du mm et aD bOOQ: Harper et al. 2005). The Euler equations are 
solved in three dimensions with a nearly isothermal polytropic 
equation of state with y = 1.1, resulting in a thermal structure 
that comes close to that obtained f rom photoionization models 
of related symbiotic binary s ystems (Nussbau mer & Vogelll 987: 
iNussbaumer & Waldedl 19931) . 

The simulation is carried out in an Eulerian frame of ref- 
erence with the stars moving within the computational domain 
which, which measures 10 15 cm a side. The computational grid 
consists of nine nested levels of refinement (Fig.[T]). The coars- 
est level consists of 50 cells cubed. From one level to the next, 
grid cells are refined by a factor of two. Each level comprises 
between 8 and 64 individual grids, the entire mesh consists of 
233 grids and 2 • 10 7 cells. The decomposed grid structure is 
exploited for parallelization under OpenMP. 

We impose free outflow at the outer domain boundary. The 
WD is represented as a spherical low pressure, low density, zero 
velocity region of radius R = 2.2 • 10 11 cm. The accretion rates 
of mass, momentum, energy, and angular momentum onto the 
WD are derived from the flow quantities entering the region. 
The rates are not sensitive to the size, pressure, and density of 
the region representing the WD. The RG is represented as a 
spherical region of 150 R©. In the rest frame of the RG star, 
which rotates and orbits with respect to the computational do- 
main, we chose a terminal velocity of the RG wind of 20 km/s. 
The absolute value of the density in the cells representing the RG 
was adjusted to achieve the desired mass loss of lO~ 7 M yr _1 . 
The wind temperature is set to 8000 K, a value between the RG 
photospheric temperature and the wind temperature more close 
to t he WD as computed by e laborate photo-ionization mod- 
els (INussbaumer & Vogell[l987h . 

At the beginning, the entire computational domain was filled 
with the RG wind. The accreting system was then relaxed over 
seven orbital periods, about the crossing time from the RG to 
the domain boundary at the RG wind speed. Mass and angular 
momentum losses out of the system were then computed through 
spherical shells around the centre of mass. 

The losses became constant for shell radii larger than 
10 14 cm. The nova explosion was then launched into the relaxed 
density structure. 
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Fig. 2. Density structure during the accretion phase. Shown is 
density (logarithmic scale, g/cm 3 ) in the orbital plane (xy-plane, 
left), and in a plane perpendicular to the orbit (yz-plane, right) 
for the entire computational domain (10 15 cm a side, top) and a 
zoom around the accreting WD (« 10 13 cm a side, bottom). Self- 
interacting high-density spirals dominate the inner region up to 
a few times the binary separation. The RG is shown in red, the 
WD in blue. 



2.2. Quiescent phase: density spirals and wind accretion 

The inter-outburst accretion phase leads to a circumstellar den- 
sity distribution that deviates substantially from the 1/r 2 density 
distribution of a single RG wind out to several system separa- 
tions. The deviations, usually neglected when interpreting RS 
Oph observations, arise from the orbital motion of the two com- 
ponents and the accretion onto the WD, and are subsequently 
transported outward in the flow. The resulting density structure, 
Fig. [2 depends critically on the ratio R = VRo/vorb, where vorb = 
2na/P, with P the orbital period. For R » 1, an Archimedian 
spiral occurs in the orbital plane, as observed i n the collid- 
ing wind Wolf Rayet binary star system WR 104 dTuthill et al.l 
1999). For RS Oph, with R < 1, the spiral pattern becomes self- 
interacting. 

In the vicinity of the WD, a high density accretion disk forms 
that is clearly visible in Fig. [2] Its diameter is roughly 1/10 of 
the system separation, and velocities are non-Keplerian. Strong 
spiral shock waves are responsible for the transport of angular 
momentum and mass, and stable accretion occurs. The disk is 
geometrically thick, the opening angle perpendicular to the or- 
bital plane (Fig. 12 lower right panel) measures roughly ±60° in 
the direction toward the RG and ±30° away from it. 

On length scales less than the system separation, radial den- 
sity profiles as seen from the WD (Fig. [3]) show density contrasts 
between the orbital plane and poleward directions of one to two 
orders of magnitude. Most prominent here, and reaching furthest 
out, is the trailing accretion wake of the WD. This density struc- 
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Fig. 3. Density profiles along radial rays originating from the 
WD. Bold colored lines denote rays in the orbital plane along the 
x- and y-coordinate axis, the RG being located in the negative x- 
direction (-x green; +x magenta; -y blue; +y red). Perpendicular, 
green bold lines denote the position of the RG. The bold black 
line is in poleward direction (+90°), thin grey lines belong to an 
inclination of 60° and positive and negative x- and y-directions, 
respectively. 



ture affects the early evolution of the nova remnant. At greater 
distances, densities in the orbital plane remain larger than toward 
the poles, but only by a factor of 2-3. This is, however, sufficient 
to cause an asymmetric evolution of the nova remnant. Note that 
although the density decreases as 1/r 2 on average, the relative 
amplitude of local variations due to the self-interacting spiral 
pattern are up to a factor of 10, especially in the orbital plane. 
We expect these local variations to vanish at even larger scales 
that are beyond our computational domain. 

The velocity of the systemic outflow is a superposition of 
the velocity of the RG wind, the stellar orbital and rotational 
velocities and hydrodynamical effects. On larger scales than 
« 5 • 10 13 cm, the flow field is essentially directed radially out- 
ward. Depending on the exact angle of an observer, any radial 
speed between 20 km/s (polewards) and « 50 km/s (mostly in 
the orbital plane) can be found in the computational data. This 
range is consistent with the different values of the wind speeds 
derived from observatio ns, e.g . 36 km/s by Iijimal (120081) and 40- 
60 km/s bv lShore et all d 19961) . 

The computed accretion rate is 10 % of the mass loss rate of 
the RG, independent of the absolute mass loss of the RG. To our 
knowledge, this is the first quantitative self-consistent prediction 
of the accretion rate of a symbiotic-like recurrent no va. A higher 

accretion rate is expected for a less massive RG (Fekel et al. 

i ^i i ii k 

2000; Dobrzycka & Kenyon 1994), since a smaller system sepa- 
ration would result. The computed value is in line with the values 
found in 3D simulations of related sy mbiotic binary star sys- 
tems, for example 6% in RW Hydrae (Dumm et al. 2000) and 
10% - 25% in Z Andromeda (M itsumoto et al.ll20Q5h 71n abso- 
lute terms, the WD in our simulation accretes lO" 8 M yr _1 . This 
value also agrees well with th e accretion rates re quired by multi- 
cycle nova evolution models (lYaron et al. 2005) for RS Oph like 
binary systems, with Mwd = 1.4M© and a recurrence period of 
22 years. Conversely, the nova models together with the above 
accretion rate of 10% of the RG mass loss quantitatively con- 
strain the mass loss of the RG to values around lO" 7 M yr _1 . 
This is yet another, completely independent, estimate of the mass 
loss rate of the RG in RS Oph. It is more at the upper limit of 
the values given by Seaquist & Taylor (1990) for RG in classical 
sy mbiotic s. A much higher estimate of 10" 5 M o yr _1 was derived 
bv lShore et al.l (Il996l) . A more recent analysis indicates the rate 
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Fig. 4. Density structure of the nova remnant. Shown are den- 
sity (logarithmic scale, g/cm3) and velocity (cm/s) in the orbital 
plane (xy-plane, left) and in a plane perpendicular to the orbit 
(yz-plane, right) at 29 hours (top panels) and 21 days (bottom 
panels) after explosion. The RG is shown in red, the WD in blue. 



could be a factor of about 10 lower (Shore 2008) and should only 
be taken as an upper limit. However, evolutionary models of sin- 
gle stars predict values greater th an lO^Mpyr" 1 for onl y a very 
short time on the RG branch (iMaeder & Mevnetlll988l) . At the 
moment, we are not able to resolve this spread in the estimates. 



2.3. Outburst: bipolar blast due to density stratification 

Launching the ejecta of a nova explosion into this complex cir- 
cumstellar environment naturally results in a strongly asymmet- 
ric shock front, visible as a high density shell in Fig |4l The 
time evolution of the nova remnant is shown in a Supplementary 
Movie. The explosion was simulated by instant injection of 
an energy and mass of 2.2 • 10 43 ergs and 2 • 10" 7 M o , re- 
spectively, at the position of the WD. The assumed explosion 
mass is on the lower end of published estimates for the ejected 
mass (ISokoloski et al.ll2006l: lO'Brien et al.lll992h . which range 
from 10" 7 M o to 10" 6 M o , but is comparable to the 2.2 • 10" 7 M o 
obtained from the self-consistently computed accretion phase. 
An isothermal equation of state with y - 1 . 1 is used, because in- 
clusion of detailed cooling as in recent ID models ( Vaytet et al. 
2007) is currently beyond reach in 3D simulations. 

The initial velocities of the modeled ejecta are ~ 3 500 km/s, 
within th e range of observed values (Bode et al. 2006; Das et al.l 
2006); lEvans et al] |2006|) . A larger extension of the shock 
front in poleward directions, where pre-o utburst densities are 
smallest, is consisten t with observations (lO'Brien et al.l 12006: 
Che sneau et al.ll2007|) . In the orbital plane, the shock front dis- 
plays a roughly circular shape for about 2/3 of its arc. Toward 
the RG and the former accretion wake, where densities are par- 
ticularly high, the shock front advances more slowly. 
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Fig. 5. Time evolution of outer shock position relative to the WD. 
Bold colored lines belong to directions within the orbital plane 
(green: toward the RG or +90°; blue: toward the trailing spiral 
of the WD or 0°; magenta: away from the RG or -90°; red: in 
direction of the orbital motion or -180°). The bold black line 
denotes the poleward direction. Thin lines belong to directions 
within the orbital plane, into the trailing accretion wake of the 
WD (50° to 60°, solid black to gray lines) and out of the system 
(-35°, dashed orange line). 



Quantitative evaluation of the shock position as a function of 
time, Fig \5\ yields shock velocities scaling with time as v^ oc t~ a , 
with 0.2 < a < 0.5 for most times and directions. Poleward 
values of a are slightly smaller after about day 10, while a = 0.5 
for directions toward the former accretion wake of the WD and 
times before day 6 after the explosion. Values for a derived from 
observational data of the 2006 nova outburst (JBode et al.l [2QQ6: 
Sok oloski et al.ll200q: iDas et al.ll2006l) are in the same range up 
to day 6, and are somewhat larger (0.45 < a < 0.79) later on. The 
present simulation suggests that the observed range of values is 
real and a consequence of different observational diagnostics 
probing different regions of the expanding shock front. A clear 
change in the velocity of the shock as a function of time is seen 
only toward the former accretion wake. 

There is no evident change in the expansion when the ac- 
cumulated mass equals the ejected mass of 2 • 10" 7 M o , which 
occurs after about eight to ten days. Observations taken much 
later will show only mixed ejecta with the nova processed mate- 
rial significantly diluted. 

2.4. Shrinking of the binary orbit 

The orbit shrinks at a rate of a/ a « -2.8 • 10" 8 per year, or 
3% per million years. This finding is in line with p revious sug- 
gestions (lHachisu et al.lll999c iJahanara et al.ll2005|) . The abso- 
lute losses of mass and angular momentum from the system, Ms 
and Js, scale linearly with the RG mass loss. For a RG mass 
loss rate of 1O" 7 M per year, the fractional losses per year are 
Ms/Ms = -2.2 • 10" 8 and J s /J s = -4.4 • 10" 8 , respectively, 
implying a dimensionless specific systemic angular momentum 
loss of (Js/Js)/(Ms/Ms) = 2. This time scale does not become 
significantly larger if the effect of the nova blast on orbital pa- 
rameters is taken into account, provided not much more mass is 
ejected than was accreted during the 22 years. 



3. Summary and conclusions 

Our 3D simulations of the quiescent phase result in a density dis- 
tribution of the circumstellar matter which is strongly stratified 
in poleward directions, leading to an expansion of the nova rem- 
nant that is roughly twice as fast in the poleward directions than 



within the orbital plane. With a reasonable RG mass loss rate 
of about 10" 7 M o yr _1 , the computed accretion rate of about 10% 
of the RG mass loss rate results in accretion rate of lO" 8 M yr _1 
fo r the WD. Ins erting this value into the multi-cycle nova models 
by lYaron et al . (2005) predicts on a qualitative level - though not 
in all quantitative details - correctly the recurrence time, ampli- 
tude, and fastness of the RS Oph nova. Our simulations further 
predict shrinking of the binary orbit at a rate of about 3% per mil- 
lion years, which would improve condi tions for enhanced mass 
transfer. According to nova models by lYaron et all (120051) . en- 
hanced mass loss should lead to shorter nova cycles and favour 
net mass gain over multiple nova cycles, a necessary condition 
for the evolution toward an SN la. On the other hand, the above 
RG mass loss suggests that the mass and thus the Roche volume 
of the RG will drastically shrink on a time scale of « 10 6 years, 
depending on the still not fixed mass of the star. The system 
then would likely enter a common envelope phase and produce 
a narrow WD-WD system. Further angular momentum losses 
by gravitational wave emission will drive the merger of the two 
WDs, an event which should be detectable with present GW de- 
tectors, e.g. Virgo and LIGO. Depending on what the result of 
the merger is, an accretion induced collapse or an SN la will end 
the life of the system RS Oph. 
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